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Abstract 

It is well known that the rotating inviscid accretion flows with adequate injection 
parameters around black holes could form shock waves close to the black holes, after 
the flow passes through the outer sonic point and can be virtually stopped by the 
centrifugal force. We examine numerically such shock waves in ID and 2D accretion 
flows, taking account of cooling and heating of the gas and radiation transport. The 
numerical results show that the shock location shifts outward compared with that 
in the adiabatic solutions and that the more rarefied ambient density leads to the 
more outward shock location. In the 2D-flow, we find an intermediate frequency 
QPO behavior of the shock location as is observed in the black hole candidate GRS 
1915+105. 

Key words: accretion, accretion disks — black hole physics — shock waves — 
radiation — hydrodynamics 



1. Introduction 

Rotating inviscid and adiabatic accretion flow around a black hole under a pseudo- 
Newtonian potentials can have two saddle type sonic points. After the inviscid flow with 
adquate injection parameters passes through the outer sonic points, the supersonic flow can be 
virtually stopped by the centrifugal force, forming a standing shock close to the black hole and 
again falling into the black hole supersonically. @Apart from the initial works in the spherical 
transonic problems of accretion and wind, it was Fukue (1987) under a full relativistic treatment 
and Chakrabarti and his collaborators (Chakrabarti 1989 , Abramowicz, Chakrabarti 1990, 
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Chakrabarti, Molteni 1993) under a pseudo-Newtonian potential who were the first to provide 
the satisfactory analytical or numerical global shock solutions for transonic, invicid, rotating 
accretion around a Schwartzschild black hole. It has been shown that these generalized accretion 
flows could be responsible for the hard and soft state transitions or the quasi-periodic oscillation 
(QPO) of the hard X-rays from the black hole candidates (Molteni et al. 1996a, Ryu 1997, 
Lanzafame et al. 1998). 

Recently, Das et al. (2001) examined analytically locations of sonic points and stand- 
ing shocks in a thin, axisymmetric, adiabatic flow around a black hole and Das (2002) and 
Das (2003) showed that, using four available pseudo-Newtonian potentials, the standing shocks 
are essential ingredients in the multi-transonic black hole accretion disks. Using the generalized 
multi-transonic accretion model. Das et al. (2003a) and Das et al. (2003b) calculated analyti- 
cally the QPO frequency of galactic black hole candidates in terms of dynamical flow variables 
and proposed a non-self-similar model of coupled accretion-outflow system in connection to 
QPO of the black hole powered galactic microquasars. Most of these analytical or numerical 
works of shocks in the rotating accretion flows around the black holes have been examined 
under the adiabatic condition, that is, the cooling and heating of gas and the radiation trans- 
port have not been taken account of. However, the radiation actually will play an important 
role in the overall flow structure, especially the shock location and the shocked temperatures. 
Accordingly it is worth-while to examine the shocks in the transonic black hole accretion disks, 
taking account of the radiation effects. 

2. Model Equations 

We examine these shock problems firstly in one-dimension and secondly in axisymmetric 
two-dimensional inviscid flow. A set of relevant equations consists of six partial differential 
equations for density, momentum, and thermal and radiation energy. These equations include 
the heating and cooling of gas and radiation transport. The radiation transport is treated in 
the gray, flux-limited diffusion approximation (Levermore, Pomraning 1981). Using cylindrical 
coordinates {r,z,ip), the basic equations for mass, momentum, gas energy, and radiation energy 
are written in the following conservative form: 




-|- div(pv) = 0, 



(1) 



d{pv) 



+ div{pvv) 




(2) 





(4) 



dpe 



+ div{pev) = —p divv — A, 



(5) 



dt 
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and 

^ + divFo + dwivEo + vPo)=A- p-^^^^v ■ Fo, (6) 
ot c 

where p is the density, v = {v,w,v^) are the three velocity components, G is the gravitational 
constant, M^, is the central mass, p is the gas pressure, e is the specific internal energy of the 
gas, Eo is the radiation energy density per unit volume, and Pq is the radiative stress tensor. It 
should be noticed that the subscript "0" denotes the value in the comoving frame and that the 
equations are correct to the first order of v/c (Kato et al. 1998). We adopt a pseudo- Newtonian 
potential (Paczyhsky, Wiita 1980) in equations (2) and (3), where Vg is the Schwartzschild radius 
given by 2GM^/c^. The force density /j^ = {fr,fz) exerted by the radiation field is given by 

fn = P^F,, (7) 

where k and a denote the absorption and scattering coefficients and Fq is the radiative fiux in 
the comoving frame. For one- dimensional form, we put z=0 in equation (2) and omit equation 
(3). 

The quantity A describes the cooling and heating of the gas, i.e., the energy exchange 
between the radiation field and the gas due to absorption and emission processes, 

A = pcK{S,-Eo), (8) 

where is the source function and c is the speed of light. For this source function, we assume 
local thermal equilibrium S^, = aT'^, where T is the gas temperature and a is the radiation 
constant. For the equation of state, the gas pressure is given by the ideal gas law, p = RcpT/ p, 
where p is the mean molecular weight and Rq is the gas constant. The temperature T is 
proportional to the specific internal energy, e, by the relation p = {'-/ — l)pe = R^pT / p, where 
7 is the specific heat ratio. To close the system of equations, we use the fiux-limited diffusion 
approximation (Levermore, Pomraning 1981) for the radiative flux: 

Fo = -— ^^gradEo, (9) 
p{K + a) 

and 

Po = Eo- TEdd, (10) 

where A and TEdd are the flux-limiter and the Eddington Tensor, respectively, for which we use 
the approximate formulas given in Kley (1989). The formulas fulflU the correct limiting condi- 
tions in the optically thick diffusion limit and the optically thin streaming limit, respectively. 

3. Model Parameters and Numerical Methods 

For the central black hole, we assume a Schwartzschild black hole with mass = IOMq. 
We try to flnd the steady solutions with shocks by solving the time-dependent equations (l)-(6). 
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Table 1. Injection Parameters 



Case 


"out 


lout 


'^out 


£out 


Rin/Rg 


Rout/rg (p 


ID 


0.0751 


0.0564 


1.875 


0.0075 


2.0 


100 


2D 


0.0938 


0.0738 


1.64 


0.005 


1.5 


30 29° 



which are numerically integrated by a finite-difference method under initial and boundary con- 
ditions. To do this, we have to determine the injection parameters such as the specific angular 
momentum Aout, the radial velocity Vont, and the sound velocity aout at an outer boundary radius 
i?out! whose parameters can lead to a shock wave close to the black hole. We search analytically 
these injection parameters through the examination of the paremeter space (e, A), where e and A 
are the total specific energy and the specific angular momentum ( e.g. Chakrabarti 1989,Molteni 
et al. 1994, Molteni et al. 1999). Typical injection parameters for ID and 2D fiows obtained 
thus are listed in table 1. Here, the velocities and distances are given in units of the speed of 
light c and the Schwartzschild radius Tg, respectively, (p is the subtended angle of the centeral 
star to the initial disk at r = -Rout, that is, tan = {h/r)ont, where h is the disk thickness. 2D 
case in the table corresponds to a super-Eddington accretion with an accretion rate M ~ 10^° 
g s~^ which is ~ 64Me for the Eddington luminosity Le = Mec^, where Me is the Eddington 
critical accretion rate, and the ambient density pout is taken to be 1.25 x 10~^ g cm~'^. 

At the outer boundary radius the injection parameters are always kept constantly. At the 
inner boundary radius -Rin, a supersonic radial velocity is specified. Adequate initial conditions 
of isothermal temperature and adequate r-dependent distribution of density are imposed on the 
flow. With these initial and boundary conditions, we perform time integration of equations (1)- 
(6) until the steady solutions are obtained. The numerical schemes used are basically the same 
as that described previously (Kley 1989,Okuda et al. 1997) but with no viscosity description. 

4. Numerical Results 

4.I. One-Dimensional Flows 

Firstly, under the injection parameters in Table 1, we numerically obtained a steady 
adiabatic shock with an ambient density pout = 3.5 x 10~^ g cm""^ at the outer boundary r = 
Rout- If a set of the injection parameters are given, the shock position in the adiabatic case is 
independent on the ambient density, because the basic equations of the adiabatic and inviscid 
flow are almost independent on the density as far as the ideal gas law is used. However, in the 
non-adiabatic cases, we obtain steady or sometimes nonsteady shock structure dependent on 
the ambient density and the shock locations are also dependent on the ambient density. Figure 
1 show the numerical results of density (a), temperature (b), radial velocity (c), and Mach 
number (d) of the velocity in the adiabatic flow (A) and non-adiabatic flows with different 
ambient densities of Pout = 3.5 X 10~5 (B), 3.5 X 10^6 (C), 7.0 x 10"^ (D) g cm'^, respectively. 
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Fig. 1. Flow structures in one dimensions with shock waves: (a) density p: red hne (A) shows the structure 
in adiabatic case with pout = 3.5 x 10~^ g cm^'^, and green (B), bhie (C), and dark blue (D) lines denote 
the structures of non-adiabtic case with pout = 3.5 x 10~^,3.5 x 10~^, and 7.0 x 10~^, (b) temperature T: 
case A is not included here, because the temperatures in case A are very hot as ^ 10^° - 10^^ K, (c) radial 
velocity v, (d) Mach number of the velocity v. 

The numerical results in the adiabatic flow agree well with the analytical ones, where the shock 
locates at r = 9.4rg. 

The shocked temperatures in the adiabatic flow are too high to be ~ 10^'^ - 10^^ K, 
while actual postshock temperatures in B, C, and D are rather low as 2 - 5 xlO^ K. The large 
temperature differences are attributed to the differences of the temperatures specifled initially 
at the outer boundary, where same injection parameters of the sound velocity aout and the 
radial velocity fout are used in the adiabatic and non-adiabatic cases. If an optically thick and 
radiation-pressure dominant accretion flow is considered in the non-adiabatic case, we have 

where Tad and T^ad are the temperatures at the outer boundary in the adiabatic and non- 
adiabatic cases, respectively. Therefore, for a given aout? we have smaller Tnad for smaller 
ambient density pout and, under the injection parameter Oout = 0.0564 in ID-flow with pout = 
3.5 X 10"^ g cm"3. Tad ~ 2 X 10^° K and Tnad ~ 1-4 x 10^ K. Tad is three orders of magnitude 
larger than Tnad- Although the temperature proflle in the non-adiabatic case is dependent 
on the ambient density, the radial velocities and their Mach numbers take the radial proflles 
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independent on the ambient density in the pre-shock region, as is found in Figure l-(c) and 
(d). 

Figure 1 also shows that the more rarefied ambient density the fiow has, the shock 
position moves more outward and the shock thickness becomes broader. The shock locations 
shift to 10.6, 12.0, and ~ 15, in cases B, C, and D, respectively. If the ambient density pout 
is too low as ~ 3.5 x 10"^ g cm~'^, the steady shock does not exist finally. At the Rankine- 
Hugoniot relation of a standing shock, the pressure balance is supported by the dominant 
radiation pressure and the ram pressure in the present non-adiabatic case. If the ambient 
density is taken to be smaller than that at a standing shock, the pressure at the upstream 
just before the shock becomes much smaller because of the lower temperature in the pre-shock 
region. Therefore, to set up a new pressure balance condition at the shock, the shock must shift 
outward as far as the same injection parameters are concerened. When the shock position shifts 
more outward, the Mach number of the pre-shock velocity decreases, that is, the shock becomes 
weaker. Generally the shock thickness is inversely proportional to the wave strength, and the 
scale factor is the radiation mean free path in our case. Therefore the shock broadens gradually 
with decreasing ambient densities. If the ambient density is too low, it may be imposibble 
to establish the Rankine-Hugoniot shock conditions in the region considered and the shock 
disappears. 

4-2. Two-Dimensional Flows 

Figures 2 and 3 show the temperature contours for 2D adiabatic and nonadiabatic ac- 
cretion fiows, respectively, under the injection parameters in Table 1 and the ambient density 
Pout = 1.25 X 10^^ g cm~^ at the outer boundary. 

There are generally two important features of a two-dimensional rotating accretion fiow 
around a black hole. One is the funnel wall, which is roughly characterized by a surface {xf,Z{) 
of vanishing effective potential 

where rf = {xf + -Zf )^^^- The other is the centrifugal barrier {xc{, z^), which is governed by the 
competition between the centrifugal force and gravitational force 

In the adiabatic fiow, the standing shock appears at r/r^ = 9.3 on the equatorial plane and is 
bent upward toward the upstream, roughly following the contour of the centrifugal barrier. The 
post-shock temperatures are as high as ~ 10^° - 10^^ K as is found in ID case. In the red and 
yellow regions between the rotating axis and the funnel wall, the temperatures are much higher 
and the densities are very low and in these regions the outfiow gas with high velocities are 
formed. These results are similar to that by the previous TVD and SPH simulations (Molteni 
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R/Rg 



Fig. 2. Temperature contours in two-dimensional adiabatic flow with shock wave. The steady shock is 
formed at r/rg 9.3 near the equatorial plane and the shock front indicated by the thick black lines 
extends obliquely to the upstream. 




R/Rg 



Fig. 3. Same as figure 2 but here the cooling and heating of the gas and the radiation transport are 
taken account of. In this figure, the shock appears at r jr^ ^ 15 near the equatorial plane but is not steady, 
and its position oscillates between 12 < rjr^ < 17. 
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R/Rg 



Fig. 4. Profiles of tfic temperature and tiie Macli number of tfie radial velocity on the equatorial plane 
for 2D adiabatic (dotted lines) and nonadiabatic (solid lines) flows. 

et al. 1996b) in the 2D adiabatic flow. 

In contrast to the adiabatic flow, when the cooling and heating of gas and the radiation 
are taken account of, a spheroidal shock is formed between the funnel wall and the equatorial 
plane and the shock oscillates quasi-periodically. Figure 3 shows the temperature contours 
at the evolutionary time t = 5.5 x lO^Vg/c in the non-adiabatic flow where the shock appears 
at Ts ~ 15rg on the equatorial plane. The overall flow is optically thick except the cone-like 
funnel region between the funnel wall and the rotational axis. In the funnel region we have 
very rarefied and optically thin flow but the temperatures are not so high compared with the 
high temperatures in the viscous accretion flows where the viscously dissipative energy heats 
up considerably (Okuda 2002). As a result, a relativistically high velocity jet is not formed 
here in spite of the high input accretion rate. The temperatures in the shocked region are 1-4 
xlO^ which are three orders of magnitude smaller than that in the adiabtic flow. For clearness 
of the shock location and the shocked temperature, in figure 4 we plot the temeprature and 
the mach number of the radial velocity on the equatorial plane. It is clear that the effects of 
radiation shifts the shock position rg ~ 9rg to ~ 15rg and the shock strength weakens. 

The shock in figure 3 oscillates and this induces the luminosity fluctuations. Figure 
5 shows the shock radii on the equatorial plane and the luminosity as a function of time in 
units of Tg/c. Although the time variation of the shock position on the equatorial plane is 
complicated, we find that it oscillates quasi-periodically with a period of ~ 2000rg/c (0.2 sec) 
and its frequency u ~ 5Hz. The spheroidal shock also oscillates in the similar way but the 
luminosity shows a more complicated variation due to the convective nature and the geometry 
of 2D-flow. 
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Fig. 5. Shock location and luminosity versus as a function of time in units of r^jc in 2D nonadiabatic 
flow. 

5. Discussion and Conclusion 

In the present paper, we have numerically examined ID and 2D inviscid transonic flows, 
taking account of the cooling and heating of gas and radiation transport. As the results, we 
flnd that the location of the centrifugally driven shock drifts more outward from the black hole 
compared with that in the adiabatic flow and that the postshock temepratures are as high as 
1-4 xlO^ K which is three orders of magnitude lower than that in the adiabatic flow. When 
a set of the injection parameters such as the speciflc angular momentum A, the radial velocity 
Wout, and the sound velocity Oout at the outer boundary are specifled, the more rarefled ambient 
density the flow has, the shock location drifts more outward. When the ambient density is 
too low, the shock wave exists no longer under the injection parameters. Depending on the 
injecteion parameters, even if the shock wave exists, it becomes sometimes unstable, that is, 
the shock location oscillates. 

The transonic accretion shocks around the black holes have been applied to 
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the hard and soft spectral states observed in the black hole candidates (Chakrabarti, 
Titarchuk 1995, Chakrabarti 1997). In the shock model, the accretion disk is decomposed 
into three distinct components : (1) an optically thick Keplerian disk on the equatorial plane 
(r > Ts); (2) a sub-Keplerian optically thin halo above the disk (r > rg); (3) a hot, optically slim 
post-shock region (r < rg ~ 5 - lOrg). The hot, dense post-shock region intercepts soft photons 
from a pre-shock region and a cold outer accretion disk and reprocesses them to form high- 
energy photons. This model seems to explain well the observed properties. In this respect, our 
results of the radiative shock may have such observational appearance in black-hole candidates. 
Molteni et al. (1996a) and Lanzafame et al. (1998) have examined shock waves in the viscous 
accretion disks around black holes and showed that when the viscosity parameter a is less than 
a critical value, standing shock waves are formed but that if the viscosity is high then the shock 
wave dissapears. Also in our ID-cases, we confirmed that the standing shock is formed for 
small a and that the shock position drifts more outward for larger viscosity parameter until it 
attains to ~ 10~^. 

The shocked accretion flows have been applied to the QPO phenomena in galactic black 
hole candidates. The galactic X-ray transient source GRS 1915-1-105 exhibits various types of 
quasi-periodic oscillations with frequencies ranging from ~ 0.001 - 10 Hz to ~ 67 Hz (Morgan 
1997). Analyzing X-ray data in GRS 1915-1-105, Rao (2000) showed that a thermal-Compton 
spectrum is responsible for the QPO generation and supported the hypothesis that the QPO is 
generated by the oscillation of the shock front in the transonic accretion flow. The QPOs are 
classified as (1) the low-frequency QPO {v\, ~ 0.001 - O.OlHz), (2) the intermediate frequency 
QPO (i^i ~ 1 - lOHz), and (3) the high frequency QPO (i/r ~ 67 Hz). The QPO-like behavior 
of the shock oscillations with i/ ~ 5 Hz in our 2D-fiow may be representative of the intermediate 
frequency QPO observed in GRS 1915-1-105. The oscillation amplitude of the shock location 
could be up to 15 % of the distance of the shock wave from the black hole. This results in a 
quasi-periodic behavior of the luminosity with an amplitude variation of a factor 2. However 
we are not able to find distinctly these frequency peaks of z/j, i^l, and i^h from the power spectra 
of the luminosity curve. 

Molteni et al. (1996a) have found the QPO behaviors in viscous accretion disks with 
shock waves around a black hole with M = IO^Mq. They obtained stable shock oscillations 
with periods of 1500 - 2000 in units of Vg/c , which depends on an index of power-law cooling 
function of gas. The oscillation period tosc of 1500 units in the case of bremsstrahlung cooling 
is comparable to 2000 units in our case. On the other hand, they interpreted the oscillation 
period as the advection timescale t^dv only when it is comparable to the cooling timescale of 
the flow and estimate it as follows: 



where w+ is the post-shock velocity. If the pre-shock velocity V- at the shock is taken to be 




(14) 
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a fraction / of the free-fall velocity, tosc = R^V^ / f where R is the compression ratio of gas 
at the shock. Using rg/rg = 17, R= 7, and / =0.7 estimated in our results, we have tosc ~ 
700 units which is smaller by a few factor than actual 2000 units. If a post-shock flow of 
constant velocity = vo/R is used, where vq is a fraction of the speed of light, tosc = Rfs/vQ. 
For Ts = 17, -R = 7, and vq = 0.066, we have tosc = 1800 units (0.17 sec) and a frequency of 
~ 6 Hz which corresponds to the shock oscillation frequency in our 2D-flow. Assuming an 
outflow from the post-shock region , Chakrabarti (1999) derived a correlation between the 
shock location and the duration of the quiescence state of GRS 1915+105 with typical QPOs 
and Chakrabarti, Manickam (2000) showed that the correlation with vq = 0.066 was well fitted 
to public-domain data from RXTE. These transonic shock wave models in rotating accretion 
flows around black holes may be promissing for the explanation of the QPOs in the black hole 
candidates. Further examinations of the models are required from a 2D-simulation point of view. 
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